FUT2_Baseline_1nd model
Structure and composition
Yiming Wang
30/05/2022
Packages
pacman::p_load(
dplyr,
tidyverse,
readxl,
vegan)
# setwd("D:PhD/01_Data/03_Vitro/02_Identification of glycan utilizing microbe") # Windows
setwd("/Volumes/Yiming_Wang/PhD/01_Data/03_Vitro/02_Identification of glycan utilizing microbe") # Mac
# setwd("/Users/godric/Desktop/Fut2 review comments/Revision") #Macbook pro
df_all <- read_excel("Glycan utilization_L6_taxa formatted.xlsx") %>%
mutate(Genotype=factor(Genotype,levels = c("WT","KO"))) %>%
mutate(Sex = factor(Sex, levels = c("Male", "Female"))) %>%
mutate(Genotype_Sex = factor(Genotype_Sex,levels = c("WT Male","WT Female","KO Male","KO Female"))) %>%
filter(Genotype == "WT") %>%
mutate(Sample = row_number()) %>%
mutate(Treatment = ifelse(Treatment == "No glycan", "mBasal", "mBasal + 2'-FL")) %>%
mutate(Treatment_Sex= paste0(Treatment,"_", Sex)) %>%
relocate(Sample, Treatment_Sex)
# Change 1 to 01
df_all$Sample <- sprintf("%02d", as.numeric(df_all$Sample))
df_all<- df_all %>%
mutate(Sample = factor(paste0("Sample", " ",df_all$Sample)))
# Subgroup by variables
#30436d WT Male
#9ea5c2 WT Female
#654e3a KO Male
#bcaba6 KO Female
#3C5488B2 WT
#7E6148B2 KO
#D98594FF Female
#94C5CCFF MaleAlpha diversity_visualization and analysis
Plot
# load package
pacman::p_load(
ggplot2,ggsci, ggthemes,
)
set.seed(123)
# Level orders
level_order_Treatment <- c("mBasal + 2'-FL", "mBasal")
# My color panel
mypal = pal_npg("nrc", alpha = 0.7)(9)
mypal## [1] "#E64B35B2" "#4DBBD5B2" "#00A087B2" "#3C5488B2" "#F39B7FB2" "#8491B4B2"
## [7] "#91D1C2B2" "#DC0000B2" "#7E6148B2"library("scales")
show_col(mypal) # #3C5488B2(WT) and #E64B35B2 (KO)
# Treatment
## FDP Diversity
Treatment_FDP <- ggplot (data=df_all, mapping = aes(y=FPD_diversity, x=factor(Treatment, level=level_order_Treatment))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Treatment),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Faith's phylogenetic diversity", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(2,10),breaks = seq(2,10,2))+
scale_colour_manual(values=c("#b86877","#70a1a7")) +
annotate("text",x=1.5,y=10,
label="p = 0.0094",hjust=0.5, size=5.5)
## Shannon Diversity
Treatment_shannon <- ggplot (data=df_all, mapping = aes(y=Shannon_diversity, x=factor(Treatment, level=level_order_Treatment))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Treatment),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Shannon diversity", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(3.5,5),breaks = seq(3.5,5,0.3))+
scale_colour_manual(values=c("#b86877","#70a1a7")) +
annotate("text",x=1.5,y=5,
label="p = 0.0020",hjust=0.5, size=5.5)
## Richness (Chao)
Treatment_chao <- ggplot (data=df_all, mapping = aes(y=Chao_richness, x=factor(Treatment, level=level_order_Treatment))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Treatment),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Chao richness", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(30,50),breaks = seq(30,50,5))+
scale_colour_manual(values=c("#b86877","#70a1a7")) +
annotate("text",x=1.5,y=50,
label="ns",hjust=0.5, size=5.5)
## Pielou evenness
Treatment_pielou <- ggplot (data=df_all, mapping = aes(y=Pielou_evenness, x=factor(Treatment, level=level_order_Treatment))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Treatment),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Pielou evenness", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(0.7,0.9),breaks = seq(0.7,0.9,0.05))+
scale_colour_manual(values=c("#b86877","#70a1a7")) +
annotate("text",x=1.5,y=0.9,
label="p < 0.0001",hjust=0.5, size=5.5)
# Print
Treatment_FDP
Treatment_shannon
Treatment_chao
Treatment_pielou
Beta diversity_results and visualisation
Datashape
# places that are not satisfied
## 1) Permanova results table
## 2) Figure A and B did not align because of legend (Change legend size, position to bottom, change P value position)
#library packages
pacman::p_load(
data.table,
vegan,
knitr,
EcolUtils
)
# Shape data
# Level orders
df_all$Treatment <- factor (df_all$Treatment, levels = c("mBasal + 2'-FL", "mBasal"))
df_all$Treatment_Sex <- factor (df_all$Treatment_Sex, levels = c("mBasal + 2'-FL_Male","mBasal + 2'-FL_Female", "mBasal_Male", "mBasal_Female"))
# subset dateframe
df_abundance <- subset(df_all, select = -c(1:15))
df_metadata <- subset(df_all, select = c(1:15))Bray-curtis
# Check the stress and choose the optimal number of dimensions
dist_bray <- vegdist(df_abundance, method = "bray")
# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)## Run 0 stress 0.01751042
## Run 1 stress 0.02006365
## Run 2 stress 0.01758811
## ... Procrustes: rmse 0.02845926 max resid 0.08098576
## Run 3 stress 0.01751045
## ... Procrustes: rmse 0.0004187265 max resid 0.001625643
## ... Similar to previous best
## Run 4 stress 0.0190382
## Run 5 stress 0.01962485
## Run 6 stress 0.01655679
## ... New best solution
## ... Procrustes: rmse 0.01742603 max resid 0.04483727
## Run 7 stress 0.01792821
## Run 8 stress 0.01704967
## ... Procrustes: rmse 0.01017602 max resid 0.03491165
## Run 9 stress 0.01655727
## ... Procrustes: rmse 0.0002291111 max resid 0.0009842971
## ... Similar to previous best
## Run 10 stress 0.01790897
## Run 11 stress 0.01790904
## Run 12 stress 0.01951685
## Run 13 stress 0.02049977
## Run 14 stress 0.01751057
## Run 15 stress 0.01775668
## Run 16 stress 0.01879263
## Run 17 stress 0.02018922
## Run 18 stress 0.01981532
## Run 19 stress 0.01792839
## Run 20 stress 0.01694019
## ... Procrustes: rmse 0.02756249 max resid 0.08516407
## *** Solution reachedNMDS$stress## [1] 0.01655679# PERMANOVA and Pair-wise PERMANOVA
set.seed(123)
## Genotype and sex are crossed variable, cage are nested in genotype and sex
adonis2(df_abundance ~ Treatment, data=df_metadata, permutations=9999, method="bray")## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = df_abundance ~ Treatment, data = df_metadata, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.85111 0.8794 131.25 1e-04 ***
## Residual 18 0.11672 0.1206
## Total 19 0.96783 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(df_abundance ~ Treatment_Sex, data=df_metadata, permutations=9999, method="bray") ## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = df_abundance ~ Treatment_Sex, data = df_metadata, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treatment_Sex 3 0.86871 0.89758 46.742 1e-04 ***
## Residual 16 0.09912 0.10242
## Total 19 0.96783 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1df_permanova_pairwise <- adonis.pair(vegdist(df_all[,16:57], method = "bray"), df_all$Treatment_Sex, nper = 9999, corr.method = "fdr")
kable(df_permanova_pairwise)| combination | SumsOfSqs | MeanSqs | F.Model | R2 | P.value | P.value.corrected |
|---|---|---|---|---|---|---|
| mBasal + 2’-FL_Male <-> mBasal + 2’-FL_Female | 0.0150608 | 0.0150608 | 1.9314491 | 0.1944781 | 0.1178 | 0.14136 |
| mBasal + 2’-FL_Male <-> mBasal_Male | 0.3907454 | 0.3907454 | 52.3473072 | 0.8674340 | 0.0103 | 0.01545 |
| mBasal + 2’-FL_Male <-> mBasal_Female | 0.3846677 | 0.3846677 | 45.2388265 | 0.8497337 | 0.0085 | 0.01545 |
| mBasal + 2’-FL_Female <-> mBasal_Male | 0.4782065 | 0.4782065 | 123.0189156 | 0.9389401 | 0.0082 | 0.01545 |
| mBasal + 2’-FL_Female <-> mBasal_Female | 0.4662045 | 0.4662045 | 94.6449896 | 0.9220615 | 0.0071 | 0.01545 |
| mBasal_Male <-> mBasal_Female | 0.0025384 | 0.0025384 | 0.5527043 | 0.0646233 | 0.7113 | 0.71130 |
# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$Treatment)
permutest(dispersion_bray)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.0014151 0.0014151 1.0824 999 0.356
## Residuals 18 0.0235310 0.0013073# Figures
#library packages
pacman::p_load(
ggthemes,
ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )
# Genotype
datascores <- as.data.frame(scores(NMDS))#$sites
scores <- cbind(as.data.frame(datascores), Treatment = df_metadata$Treatment)
centroids <- aggregate(cbind(NMDS1, NMDS2) ~ Treatment, data = scores, FUN = mean)
seg <- merge(scores, setNames(centroids, c('Treatment','oNMDS1','oNMDS2')),
by = c('Treatment'), sort = FALSE)
composition_Treatment <- ggplot(scores, aes(x = NMDS1, y = NMDS2, color = Treatment)) +
geom_segment(data = seg,
mapping = aes(xend = oNMDS1, yend = oNMDS2),
size=0.25) +
geom_point(data = centroids, size = 4) +
geom_point()+
coord_fixed()+
theme_bw() +
labs(title="")+
theme(legend.position = "bottom",
legend.title= element_text(size=10, face="plain"),
legend.text = element_text(size=14),
legend.key = element_rect(fill = "transparent", colour = "black"),
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 9), size=14),
axis.title.y=element_text(margin = margin(r = 5), size=14),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_x_continuous(limits = c(-0.5,0.5),breaks = seq(-0.5,0.5,0.25))+
scale_y_continuous(limits = c(-0.25,0.25),breaks = seq(-0.25,0.25,0.1))+
scale_color_manual(values=c("#70a1a7", "#b86877"))+
scale_fill_manual(values=c("#70a1a7", "#b86877"))+
annotate("text",x=-0.5,y=-0.25, #0.4 if legend position is right
label="Stress = 0.017",hjust=0, size=4.5)+ #2.2 if legend position is right
annotate("text",x=-0,y=0.25,
label=expression(italic("P") < "0.0001"),hjust=0.5, size=4.5)
composition_Treatment
ANOISM: Variance between groups and within groups
# Bray-curtis
dist_bray <- vegdist(df_abundance, method = "bray")
## Genotype
anosim.Treatment_bray <- with(df_metadata, anosim(dist_bray, Treatment))
summary(anosim.Treatment_bray)##
## Call:
## anosim(x = dist_bray, grouping = Treatment)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.998
## Significance: 0.001
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0971 0.1675 0.2408 0.3287
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 87 115.75 140.5 165.25 190 100
## mBasal + 2'-FL 1 29.00 56.0 79.00 93 45
## mBasal 3 21.00 40.0 57.00 77 45plot(anosim.Treatment_bray)